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Abstract. In this paper, we study robust estimators of the memory parameter d of a (possi- 
bly) non stationary Gaussian time series with generalized spectral density /. This generalized 
spectral density is characterized by the memory parameter d and by a function /* which 
specifies the short-range dependence structure of the process. Our setting is semi-parametric 
since both /* and d are unknown and d is the only parameter of interest. The memory pa- 
rameter d is estimated by regressing the logarithm of the estimated variance of the wavelet 
coefficients at different scales. The two estimators of d that we consider are based on robust 
estimators of the variance of the wavelet coefficients, namely the square of the scale estimator 
proposed by [27] and the median of the square of the wavelet coefficients. We establish a 
Central Limit Theorem for these robust estimators as well as for the estimator of d based on 
the classical estimator of the variance proposed by [19]. Some Monte-Carlo experiments are 
presented to illustrate our claims and compare the performance of the different estimators. 
The properties of the three estimators are also compared on the Nile River data and the In- 
ternet traffic packet counts data. The theoretical results and the empirical evidence strongly 
suggest using the robust estimators as an alternative to estimate the memory parameter d 
of Gaussian time series. 



1. Introduction 

Long-range dependent processes are characterized by hyperbolically slowly decaying corre- 
lations or by a spectral density exhibiting a fractional pole at zero frequency. During the last 
decades, long-range dependence (and the closely related self-similarity phenomena) has been 
observed in many different fields, including financial econometrics, hydrology or analysis of 
Internet traffic. In most of these applications, however, the presence of atypical observations 
is quite common. These outliers might be due to gross errors in the observations but also to 
unmodeled disturbances; see for example [31] and [30] for possible explanations of the presence 
of outliers in Internet traffic analysis. It is well-known that even a few atypical observations 
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can severely affect estimators, leading to incorrect conclusions. Hence, defining robust esti- 
mators of the memory parameter which are less sensitive to the presence of additive outliers 
is a challenging practical problem. 

In this paper, we consider the class of fractional processes, denoted M(d) defined as follows. 
Let X = {Xk}k&z be a real-valued Gaussian process, not necessarily stationary and denote 
by A.X the first order difference of X, defined by [AX] n = X n — X n _i, n £ Z. Define, 
for an integer K > 1, the ET-th order difference recursively as follows : A A = A o A. K ~ l . 
Let /* be a bounded non-negative symmetric function which is bounded away from zero in a 
neighborhood of the origin. Following [20], we say that X is an M(d) process if for any integer 
K > d — 1/2, A. K X is stationary with spectral density function 

/ A K X (A) = |l-e- iA | 2 ^- d )r(A), Ae(-7r,7r). (1) 

Observe that /a^x(^) i n (1) i s integrable since —(K — d) < 1/2. When d > 1/2, the process 
is not stationary. One can nevertheless associate to X the function 

/(A) = |l-e- iA |- 2d r(A), (2) 

which is called a generalized spectral density function. In the sequel, we assume that /* € 
%{f3,L) with < f3 < 2 and L > where H(f3,L) denotes the set of non-negative and 
symmetric functions g satisfying, for all A £ (— ir,ir), 

\g(X)-g(0)\<Lg(0)\Xf . (3) 

Our setting is semi-parametric in that both d and /* in (2) are unknown. Here, /* can be 
seen as a nuisance parameter whereas d is the parameter of interest. This assumption on 
/* is typical in the semi-parametric estimation setting; see for instance [25] and [21] and the 
references therein. 

Different approaches have been proposed for building robust estimators of the memory 
parameter for M(d) processes in the semi-parametric setting outlined above. [31] have pro- 
posed a robustified wavelet based-regression estimator developed by [1]; the robustification 
is achieved by replacing the estimation of the wavelet coefficients variance at different scales 
by the median of the square of the wavelet coefficients. Another technique to robustify the 
wavelet regression technique has been outlined in [23] which consists in regressing the loga- 
rithm of the square of the wavelet coefficients at different scales. [18] proposed a robustified 
version of the log-periodogram regression estimator introduced in [14]. The method replaces 
the log-periodogram of the observation by a robust estimator of the spectral density in the 
neighborhood of the zero frequency, obtained as the discrete Fourier transform of a robust 
autocovariance estimator defined in [17]; the procedure is appealing and has been found to 
work well but also lacks theoretical support in the semi-parametric context (note however that 
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the consistency and the asymptotic normality of the robust estimator of the covariance have 
been discussed in [16]). 

In the related context of the estimation of the fractal dimension of locally self-similar Gauss- 
ian processes [10] has proposed a robust estimator of the Hurst coefficient; instead of using 
the variance of the generalized discrete variations of the process (which are closely related to 
the wavelet coefficients, despite the facts that the motivations are quite different), this author 
proposes to use the empirical quantiles and the trimmed- means. The consistency and asymp- 
totic normality of this estimator is established for a class of locally self-similar processes, using 
a Bahadur- type representation of the sample quantile; see also [9]. [28] proposes to replace 
the classical regression of the wavelet coefficients by a robust regression approach, based on 
Huberized M-estimators. 

The two robust estimators of d that we propose consist in regressing the logarithm of robust 
variance estimators of the wavelet coefficients of the process X on a range of scales. We use 
as robust variance estimators the square of the scale estimator proposed by [27] and the 
square of the mean absolute deviation (MAD). These estimators are a robust alternative to 
the estimator of d proposed by [19] which uses the same method but with the classical variance 
estimator. Here, we derive a Central Limit Theorem (CLT) for the two robust estimators of 
d and, by the way, we give another methodology for obtaining a Central Limit Theorem for 
the estimator of d proposed by [19]. In this paper, we have also established new results on the 
empirical process of array of stationary Gaussian processes by extending [3, Theorem 4 ] and 
the Theorem of [11] to arrays of stationary Gaussian processes. These new results were very 
helpful in establishing the CLT for the three estimators of d that we propose. 

The paper is organized as follows. In Section 2, we introduce the wavelet setting and 
define the wavelet based regression estimators of d. Section 3 is dedicated to the asymptotic 
properties of the robust estimators of d. In this section, we derive asymptotic expansions of 
the wavelet spectrum estimators and provide a CLT for the estimators of d. In Section 4, some 
Monte-Carlo experiments are presented in order to support our theoretical claims. The Nile 
River data and two Internet traffic packet counts datasets collected from the University of 
North Carolina, Chapel are studied as an application in Section 5. Sections 6 and 7 detail the 
proofs of the theoretical results stated in Section 3. 
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2. Definition of the wavelet-based regression estimators of the memory 

PARAMETER d. 

2.1. The wavelet setting. The wavelet setting involves two functions cp and tp in L 2 (M) and 
their Fourier transforms 

0(£) d ^ r 0(t) e -^ dt and = ^(t)e-'^ dt . (4) 

J —oo J —oo 

Assume the following: 

(W-l) 4> and tp are compactly-supported, integrable, and (p(0) = <p{t) dt = 1 and 

n v^(t)dt=i. 

(W-2) There exists a > 1 such that sup ?eIR (1 + < oo. 

(W-3) The function tp has M vanishing moments, i.e. t m ip(t)dt = for all m = 

0,...,M- 1. 

(W-4) The function X^fcez k m (p(- — k) is a polynomial of degree m for all m = 0, . . . , M — 1. 

Condition (W-2) ensures that the Fourier transform ^ decreases quickly to zero. Condi- 
tion (W-3) ensures that tp oscillates and that its scalar product with continuous-time polyno- 
mials up to degree M — 1 vanishes. It is equivalent to asserting that the first M —1 derivatives 
of tp vanish at the origin and hence 

\ip(X)\ = 0(\X\ M ), asA^O. (5) 

Daubechies wavelets (with M > 2) and the Coiflets satisfy these conditions, see [19]. Viewing 
the wavelet ip(t) as a basic template, define the family {ipj t k,j & Z, k £ Z} of translated and 
dilated functions 

rp jtk (t) = 2~ j/2 ip(2~ j t - k), j G Z, A; G Z . (6) 

Positive values of k translate tp to the right, negative values to the left. The scale index j 
dilates ip so that large values of j correspond to coarse scales and hence to low frequencies. 
We suppose throughout the paper that 

(l + p)/2-a<d<M. (7) 

We now describe how the wavelet coefficients are defined in discrete time, that is for a real- 
valued sequence {xfc, k 6 Z} and for a finite sample {x k , k = 1, . . . ,n}. Using the scaling 
function (p, we first interpolate these discrete values to construct the following continuous-time 
functions 

n 

x n (t) A =^2x k 4){t-k) and x(t) d = ^ x k (p(t - k), teR. (8) 

k=l k&L 
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Without loss of generality we may suppose that the support of the scaling function (f> is 
included in [— T,0] for some integer T > 1. Then 

x n (i) = x(i) for all * 6 [0, n - T + 1] . 

We may also suppose that the support of the wavelet function ip is included in [0, T]. With 
these conventions, the support of tpj ^ is included in the interval [2 J k, 2 J (k + T)]. The wavelet 
coefficient Wj >k at scale j > and location k E Z is formally defined as the scalar product in 



L 2 (R) of the function 1 1— > x(t) and the wavelet t i— > ijjj k(t): 



W jtk = / x(t)^- fc (t) dt = / x n (t)^, fc (t) dt, j > 0, A; € Z , (9) 

J —oo J —oo 

when [2 J '/c, 2 J 'fc + T] C [0, n — T + 1], that is, for all (j, fe) G X n , where 

X n d J: f {(j,A;) : j > 0,0 < fc < nj - 1} with nj = [2~ j (n - T + 1) - T + 1] . (10) 

If A M X is stationary, then from [20, Eq (17)] the process {Wj t k}keZ of wavelet coefficients 
at scale j > is stationary but the two-dimensional process {[W^fc, Wj^k] T }keZ of wavelet 
coefficients at scales j and j', with j > j', is not stationary. Here T denotes the transposition. 
This is why we consider instead the stationary between-scale process 

{\Wj,k,VV j ,k(j-j') T ] T }k&, (11) 
where ~Wj t k(j — j') is defined as follows: 

def r 1 T 

W j,kU - f) = { W ji,2J-j'ki W j',23-i'k+V • • • ' W j',2J-j'k+2J~j'-l ■ 

For all j,j' > 1, the covariance function of the between scale process is given by 

Cov(W A k >(j - j'),Wj, k ) = e iA ( fc " fc ') Djj_j/(A; /) dX , (12) 

J — IT 

where Dj (A; /) stands for the cross-spectral density function of this process. For further 
details, we refer the reader to [20, Corollary 1]. The case j = j' corresponds to the spectral 
density function of the within-scale process {Wj^keZ- 

In the sequel, we shall use that the within- and between-scale spectral densities Dj-j_j/(A; d) 
of the process X with memory parameter d G R can be approximated by the corresponding 
spectral density of the generalized fractional Brownian motion defined, for d £ R and 
u £ N, by 

D OC)U (A; d) = [D W u (A; d), . . . , D^" 1 ) (A; d) 



| A + 2/^|~ M e u (A + 2/tt) ^(A + 2Z7r)^(2 _ "(A + 2Zvr)) , (13) 
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where, 

e u (0 d ^ f 2^/ 2 [1, e" i2 ^, . . . , e-^- 1 ) 2 -^, fG K. 
For further details, see [19, p. 307]. 

2.2. Definition of the robust estimators of d. Let us now define robust estimators of 
the memory parameter d of the M(d) process X from the observations X±, . . . ,X n . These 
estimators are derived from the [1] construction, and consists in regressing estimators of the 
scale spectrum 

a] = Var(W,- >0 ) (14) 

with respect to the scale index j. More precisely, if cr 2 is an estimator of <t 2 based on 
Wjfi-.rij-i = (Wjfi, . . . , Wj in -— i) then an estimator of the memory parameter d is obtained by 
regressing log(<7 2 ) for a finite number of scale indices j £ {Jo, . . . , Jo+i} where Jo = Jo(n) > 
is the lower scale and 1 + 1 > 2 is the number of scales in the regression. The regression esti- 
mator can be expressed formally as 

j +e 

d n (J ,w) = ^2 wj-, Jo log (a j) , (15) 
j=Jo 

where the vector w = f [wq, . . . , wg] T of weights satisfies ^i=o w i = ® an< ^ ^ l°g(2) Y2i=o = ^-j 
see [1] and [20]. For Jq > 1 and i > 1, one may choose for example w corresponding to the 
least squares regression matrix, defined by w = DB(B T DB)~ 1 h where 



, def 

b = 



(21og(2)) 



-i 



fi drf 



T 

1 1 ... 1 

1 ... I 



is the design matrix and D is an arbitrary positive definite matrix. The best choice of D 
depends on the memory parameter d. However a good approximation of this optimal matrix 
D is the diagonal matrix with diagonal entries Di^ = 2~ l , i = 0...,£; see [13] and the 
references therein. We will use this choice of the design matrix in the numerical experiments. 
A heuristic justification for this choice is that by [19, Eq. (28)], 



2 

<7„- ~ 



C2 2jd , as j -»• oo , (16) 



where C is a positive constant. 

In the sequel, we shall consider three different estimators of d based on three different 
estimators of the scale spectrum <r| with respect to the scale index j which are defined below. 
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2.2.1. Classical scale estimator. This estimator has been considered in the original contribu- 
tion of [1] and consists in estimating the scale spectrum crj with respect to the scale index j 
by the empirical variance 

rij 

-c Lj = ^£^, (17) 

7 ■ i 
J 1=1 

where for any j, rij denotes the number of available wavelet coefficients at scale index j defined 
in (10). 

2.2.2. Median absolute deviation. This estimator is well-known to be a robust estimator of the 
scale and as mentioned by [27] it has several appealing properties: it is easy to compute and 
has the best possible breakdown point (50%). Since the wavelet coefficients Wjj are centered 
Gaussian observations, the square of the median absolute deviation of Wj^o-.nj-i is defined by 

°MADj = (™(®) ™ ed J W j,i\) > ( 18 ) 
\ 0<i<rij — 1 J 

where $ denotes the c.d.f of a standard Gaussian random variable and 

m($) = l/$ _1 (3/4) = 1.4826 . (19) 

The use of the median estimator to estimate the scalogram has been suggested to estimate 
the memory parameter in [29] (see also [24, p. 420]). A closely related technique is considered 
in [9] and [10] to estimate the Hurst coefficient of locally self-similar Gaussian processes. Note 
that the use of the median of the squared wavelet coefficients has been advocated to estimate 
the variance at a given scale in wavelet denoising applications; this technique is mentioned 
in [12] to estimate the scalogram of the noise in the i.i.d. context; [15] proposed to use this 
method in the long-range dependent context; the use of these estimators has not been however 
rigorously justified. 

2.2.3. The Croux and Rousseeuw estimator. This estimator is another robust scale estimator 
introduced in [27]. Its asymptotic properties in several dependence contexts have been further 
studied in [16] and the square of this estimator is defined by 

£crj = (c(*){|W,-i - W jjk \; 0<i,k< nj - l} (fen .)) 2 , (20) 

where c($) = 2.21914 and k nj = [nf/4j- That is, up to the multiplicative constant c(<&), 
&CR,j is the k nj th order statistics of the distances \Wj t i — Wj^\ between all the pairs of 
observations. 
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3. Asymptotic properties of the robust estimators of d 



3.1. Properties of the scale spectrum estimators. The following proposition gives an 
asymptotic expansion for &chj> ^MAD j ano ^ &CRj defined in (17), (18) and (20), respectively. 
These asymptotic expansions are used for deriving Central Limit Theorems for the different 
estimators of d. 

Proposition 1. Assume that X is a Gaussian M(d) process with generalized spectral density 
function defined in (2) such that f* G H(/3, L) for some L > and < (3 < 2. Assume 
that (W-l)-(W-4) hold with d, a and M satisfying (7). Let Wj^ be the wavelet coefficients 
associated to X defined by (9). J/nH Jo(n) is an integer valued sequence satisfying Jo(n) — > 
oo and n2" J °^ — > oo, as n — >■ oo ; then a 2 a defined in (17), (18) and (20), satisfies the 
following asymptotic expansion, as n —> oo, for any given i > 1 



max 

Mn)<j<J (n)+£ 



n j@lj ~ of) 



2<r? ' fWu 



o P (l) 



(21) 



where * denotes CL, CR and MAD 7 a 2 is defined in (14) and IF is given by 



IF(x,CL,#) = -ff 2 (x), 

IF (,, CR, •) = c(t) (I/- 4 - + 1/Cm) + * {X - ^ » 



IF(x,MAD,$) 



-m(<&) 



(l{x<l/m($)} _ 3/4) - (l{ g <-i/ m ($)} - 1/4) 
2<p(l/m($)) 



(22) 
(23) 

(24) 



where (p denotes the p.d.f of the standard Gaussian random variable, m($) and c($) being 
defined in (19) and (20), respectively and H2(x) = x 2 — 1 is the second Hermite polynomial. 

The proof is postponed to Section 6. 

We deduce from Proposition 1 and Theorem 6 given and proved in Section 6 the following 
multivariate Central Limit Theorem for the wavelet coefficient scales. 



Theorem 2. Under the assumptions of Proposition 1, (a 2 Jq , . . . , a 2 j Q+i ) T , where a 2 ■ 
defined in (17), (18) and (20), satisfies the following multivariate Central Limit Theorem 



IS 



Vn2- J o2- 2Jod 


( 






2 




AaT(0,U*(cZ)) 




\ 


_ a *,J +e _ 




2 

_ a *,Jo+l _ 


) 





(25) 
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where 

-2 

p\K(d)P 



U*,iAd) = 4(/* (0)) 2 £ ,^; } _ 2 2 d ^^ j2d (2- P )iM + iM 



p>2 

x E 2 E 1 ( f D S|^|( A ^) eiATdA ) P . < < t ■ ( 26 ) 

in (26), K((i) = f J K |Cl _2rf |V ; (0|dC) (•; d) is the cross-spectral density defined in (13), 

c p (IF*) = ~E\]F(X,*,$>)Hp(X)], where H p is the pth Hermite polynomial and IF(-,*,$) is 
defined in (22), (23) and (24). 

The proof of Theorem 2 is postponed to Section 6. 

Remark 1. Since for * = CL, IF(-) = H2O/2, Theorem 2 gives an alternative proof to [19, 
Theorem 2] of the limiting covariance matrix of (<7 2 L j , . . . , Uq L j +i) T which is given, for 
0<i,j<£, by 

V C L,iAd) = 4tt(/*(0)) 2 2 4 ^)+^ r {Tt^^d^dX . 

J —TT 

Thus, for * = CR and * = MAD, we deduce the following 

UcL,i,i(d) 1/2 . , 

U. lM (d) - M[lFl(Z)] ' 1 ; 

where Z is a standard Gaussian random variable. With Lemma 8, we deduce from the inequal- 
ity (27) that the asymptotic relative efficiency of cr 2 • is larger than 36. 76% when * = MAD 
and larger than 82.27% when * = CR. 

3.2. CLT for the robust wavelet-based regression estimator. Based on the results 
obtained in the previous section, we derive a Central Limit Theorem for the robust wavelet- 
based regression estimators of d defined by 

J +£ 

d* >n (Jo,w) = E w i~ J o lo & ' ( 28 ) 

3 = Jo 

where a 2 ■ are given for * = CL, MAD and CR by (17), (18) and (20), respectively. 
Theorem 3. Under the same assumptions as in Proposition 1 and if 

n2 -(l+2/3)J (n) _^ q ^ as n ^ (29) 

then, d* in (Jo,w) satisfies the following Central Limit Theorem: 

\/n2~ J o(n) (d %n (J ,w) - d\ A M (0, w r V*(d)w) , (30) 
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where V*(d) is the (1 + i) x (1 + £) matrix defined by 

v *« («o = E ^IS 2Hl " Jl+iAJ E 21 E 1 ( f d > iATdA ) P • * ^' ^ ■ 

p>2 P W r G Z r=0 • y " 7r 

(31) 

In (31), K(d) = J K |Cl~ 2d |V , (C)|dC) d) is i/ie cross-spectral density defined in (13), 

Cp (IF*) = E[IF(X, *, $)-Hp(X)], where H p is the pth Hermite polynomial and IF(-,*,<3?) is 
defined in (22), (23) and (24). 

The proof of Theorem 3 is a straightforward consequence of [19, Proposition 3] and Theo- 
rem 2 and is thus not detailed here. 

Remark 2. Since it is difficult to provide a theoretical lower bound for the asymptotic relative 
efficiency (ARE) o/ d* iri (Jn, w) defined by 

ARE*(d) = w T V CL (d)w/w T V # ((i)w , (32) 

where * = CR or MAD, we propose to compute this quantity empirically. We know from The- 
orem 3 that the expression of the limiting covariance matrix V *,i,j(d) is valid for all Gaussian 
M(d) processes satisfying the assumptions given in Proposition 1, thus it is enough to compute 
ARE*(d) in the particular case of a Gaussian ARFIMA(0,d,0) process (Xt). Such a process 
is defined by 

X t = (I - BY d Z t = Y pl+A- Z^ , (33) 

v > r(j + 1)r{d) j k J 

where {Z t } are i.i.d M(0, 1). We propose to evaluate ARE*(d) when d belongs to [—0.8; 3]. 
With such a choice of d, both stationary and non- stationary processes are considered. The 
empirical values of ARE*(d) are given in Table 1. The results were obtained from the obser- 
vations Xi, . . . , X n where n = 2 12 and 1000 independent replications. We used Daubechies 
wavelets with M = 2 vanishing moments when d < 2 and M = 4 when d > 2 which ensures 
that condition (7) is satisfied. The smallest scale is chosen to be Jq = 3 and Jq + £ = 8. 



d 


-0.8 


-0.4 


-0.2 





0.2 


0.6 


0.8 


1 


1.2 


1.6 


2 


2.2 


2.6 


3 


ARE C n{d) 
AREmao {d) 


0.72 
0.48 


0.67 
0.39 


0.63 
0.38 


0.65 
0.36 


0.70 
0.43 


0.63 
0.39 


0.70 
0.44 


0.75 
0.47 


0.76 
0.45 


0.75 
0.50 


0.79 
0.48 


0.74 
0.5 


0.77 
0.49 


0.74 
0.49 



TABLE 1. Asymptotic relative efficiency of rf„,CR and d n> MAD with respect to d n ,cL- 



From Table 1, we can see that d Ut cR is more efficient than d nj MAD an d that its asymptotic 
relative efficiency AREqr ranges from 0.63 to 0.79. These results indicate empirically that 
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the the loss of efficiency of the robust estimator d nt QR is moderate and makes it an attractive 
robust procedure to the non-robust estimator d n; cL- 

4. Numerical experiments 

In this section the robustness properties of the different estimators of d, namely dch,n{Jo, w )j 
^CR,n(^0) w ) an d dyiAD,n(Jo, w ) 5 that are defined in Section 2.2 are investigated using Monte 
Carlo experiments. In the sequel, the memory parameter d is estimated from n = 2 12 obser- 
vations of a Gaussian ARFIMA(0,d,0) process defined in (33) when d—0.2 and 1.2 eventually 
corrupted by additive outliers. We use the Daubechies wavelets with M = 2 vanishing mo- 
ments which ensures that condition (7) is satisfied. 

Let us first explain how to choose the parameters Jo an d Jo + i- With n = 2 12 , the maximal 
available scale is equal to 10. Choosing Jo too small may introduce a bias in the estimation 
of d by Theorem 3. However, at coarse scales (large values of Jo), the number of observations 
may be too small and thus choosing Jo too large may yield a large variance. Since at scales 
j = 9 and j = 10, we have respectively 5 and 1 observations, we chose Jo + £ = 8. For the 
choice of Jo, we proposed to use the empirical rule illustrated in Figure 1. In this figure, we 
display the estimates d nt ch, dn,CR an d dn,MAD of the memory parameter d as well as their 
respective 95% confidence intervals from Jo = 1 to Jo = 7 with Jo + t = 8. We propose to 
choose Jo = 3 in both cases (d = 0.2 and d = 1.2) since the successive confidence intervals 
starting from Jo = 3 to Jo = 7 are such that the smallest one is included in the largest one. 
We shall take Jq = 3 in the sequel. 




0246 0246 

FIGURE 1. Confidence intervals of the estimates d„,CL, d nj cR and d n ,MAD of an 
ARFIMA(0, d, 0) process with d = 0.2 (left) and d = 1.2 (right) for Jo = 1,...,8 and 
Jo+i = 9. For each Jo, are displayed confidence interval associated to d„,cL (red), 
intervald nj cR (gren) and d n ,MAD (blue), respectively. 
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In the left panels of Figures 2 and 3 the empirical distribution of Vn2~ J °(d*^ n — d) are 
displayed when * = CL, MAD and CR for the ARFIMA(0,d,0) model with d = 0.2 (Figure 2) 
and d = 1.2 (Figure 3), respectively. They were computed using 5000 replications; their 
shapes are close to the Gaussian density (the standard deviations are of course different). 
In the right panels of Figures 2 and 3, the empirical distribution of V n2~ J °(d* tn — d) are 
displayed when outliers are present. We introduce 1% of additive outliers in the observations; 
these outliers are obtained by choosing uniformly at random a time index and by adding to 
the selected observation 5 times the standard error of the raw observations. The empirical 
distribution of V n2~ J °(dcL,n — d) is clearly located far away from zero especially in the non 
stationary ARFIMA(0, 1.2, 0) model. One can also observe the considerable increase in the 
variance of the classical estimator. In sharp contrast, the distribution of the robust estimators 
Vn2~J° (dyiAD,n ~ d) and V n2~ J °(dcR, n — d) stays symmetric and the variance stays constant. 




-4 -2 2 4 -4 -2 2 



FIGURE 2. Empirical densities of the quantities \fn2 ^o(d*, n — d), with * = CL (solid 
line), * = CR (dashed line) and * = MAD (dotted line) of the ARFIMA(0,0.2,0) model 
without outliers (left) and with 1% of outliers (right). 

5. Application to real Data 

In this section, we compare the performance of the different estimators of the long memory 
parameter d introduced in Section 2.2 on two different real datasets. 

5.1. Nile River data. The Nile River dataset is a well-known time series, which has been 
extensively analyzed; see [5, Section 1.4, p. 20]. The data consists of yearly minimal water 
levels of the Nile river measured at the Roda gauge, near Cairo, for the years 622-1284 AD 
and contains 663 observations; The units for the data as presented by [5] are centimeters. 
The empirical mean and the standard deviation of the data are equal to 1148 and 89.05, 
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FIGURE 3. Empirical densities of the quantities V n2~ J o (d*, n — of), with * = CL (solid 
line), * = CR (dashed line) and * = MAD (dotted line) of the ARFIMA(0, 1.2,0) model 
without outliers (left) and with 1% of outliers (right). 

respectively. The question has been raised as to whether the Nile time series contains outliers; 
see for example [4], [25], [8] and [18]. The test procedure developed by [8] suggests the presence 
of outliers at 646 AD (p- value 0.0308) and at 809 (p- value 0.0007). Another possible outliers 
is at 878 AD. Since the number of observations is small, in the estimation of d, we took Jq = 1 
and Jo + I = 6. With this choice, we observe a significant difference between the classical 
estimators d ny cl = 0.28 (with 95% confidence interval [0.23, 0.32]) and the robust estimators 
d n>C R = 0.408 (with 95% confidence interval [0.34, 0.46]) and cf„ iM AD = 0.414 (with 95% 
confidence interval [0.34, 0.49]). Thus, to better understand the influence of outliers on the 
estimated memory parameter in practical situations, a new dataset with artificial outliers was 
generated. Here, we replaced the presumed outliers of [8] by the value of the observation plus 
10 times the standard deviation. The new memory parameter estimators are d n ,CL = 0.12, 
d n ,CR = 0.4 and d n ,MAD = 0.392. As was expected, the values of the robust estimators 
remained stable. However, the classical estimator of d was significantly affected. A robust 
estimate of d for the Nile data is also given in [2] and in [18]. The authors found 0.412 and 
0.416, respectively. These values are very close to d n QR = 0.408 and d nj MAD = 0.414. 

5.2. Internet traffic packet counts data. In this section, two Internet traffic packet counts 
datasets collected at the University of North Carolina, Chapel (UNC) are analyzed. These 
datasets are available from the website http:/ /netlab. cs.unc.edu/public/old_research/net_lrd/ 
These datasets have been studied by [23]. 

Figure 4 (left) displays a packet count time series measured at the link of UNC on April 
13, Saturday, from 7:30 p.m. to 9:30 p.m., 2002 (Satl930). Figure 4 (right) displays the same 
type of time series but on April 11, a Thursday, from 1 p.m. to 3 p.m., 2002 (Thul300). These 
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packet counts were measured every 1 millisecond but, for a better display, we aggregated them 
at 1 second. 




2000 4000 6000 2000 4000 6000 



Sat1930 Thu1300 

FIGURE 4. Packet counts of aggregated traffic every 1 second. 

The maximal available scale for the two datasets is 20. Since we have less than 4 observations 
at this scale, we set the coarse scale Jo + t = 19 and vary the finest scale Jo from 1 to 17. 
The values of the three estimators of d are stored in Table 2 for Jo = 1 to 14 as well as the 
standard errors of V n2~ J °{d nt , — d) for the two datasets: Thul300 and Satl930. 



Jo 


1 


2 


3 


4 


5 


6 


7 


8 


9 


10 


11 


12 


13 


14 


Thul300 


d n ,CL 


0.08 


0.09 


0.11 


0.15 


0.19 


0.25 


0.31 


0.39 


0.43 


0.47 


0.51 


0.49 


0.44 


0.41 


SE CL 


(0.52) 


(0.56) 


(0.51) 


(0.52) 


(0.57) 


(0.52) 


(0.56) 


(1.45) 


(0.74) 


(0.76) 


(0.87) 


(0.91) 


( i.io) 


(1.21) 


dn,CR 


0.08 


0.07 


0.07 


0.09 


0.13 


0.19 


0.28 


0.34 


0.37 


0.40 


0.42 


0.43 


0.48 


0.45 


SE C r 


(0.55) 


(0.58) 


(0.61) 


(0.63) 


(0.59) 


(0.6) 


(0.67) 


(1.42) 


(0.82) 


(0.88) 


(0.97) 


(1.08) 


(1.18) 


(1.23) 


4.MAD 


0.08 


0.08 


0.07 


0.09 


0.13 


0.19 


0.27 


0.33 


0.38 


0.40 


0.43 


0.43 


0.5 


0.48 


SEmad 


(0.74) 


(0.87) 


(0.78) 


(0.83) 


(0.86) 


(0.84) 


(0.91) 


(1.49) 


(0.98) 


(1.04) 


(1.07) 


(1.15) 


(1.18) 


(1.2) 


Satl930 


d n ,CL 


0.05 


0.06 


0.08 


0.11 


0.14 


0.17 


0.23 


0.28 


0.33 


0.36 


0.37 


0.39 


0.42 


0.42 


SE CL 


(0.41) 


(0.47) 


(0.43) 


(0.48) 


(0.47) 


(0.48) 


(0.46) 


(0.89) 


(0.54) 


(0.61) 


(0.70) 


(0.80) 


(1.11) 


(1.24) 


dn,CR 


0.06 


0.06 


0.06 


0.09 


0.12 


0.16 


0.23 


0.3 


0.34 


0.38 


0.4 


0.42 


0.44 


0.42 


SE CR 


(0.51) 


(0.47) 


(0.54) 


(0.48) 


(0.48) 


(0.53) 


(0.56) 


(0.90) 


(0.81) 


(0.70) 


(0.88) 


(0.96) 


(1.21) 


(1.26) 




0.06 


0.06 


0.07 


0.09 


0.11 


0.16 


0.23 


0.29 


0.33 


0.38 


0.4 


0.43 


0.45 


0.4 


SEmad 


(0.59) 


(0.77) 


(0.72) 


(0.81) 


(0.70) 


(0.89) 


(0.82) 


(0.64) 


(1.13) 


(0.99) 


(1.10) 


(1.34) 


(1.49) 


(1.38) 



TABLE 2. Estimators of d with J = 1 to J = 14 and J +£ = 19 obtained from Thul300 



and Satl930. Here SE denotes the standard error of v n2~ J o (d„,» — d). 



In Figure 5, we display the estimates d nj cL> d nt cR and d nj MAD of the memory parameter d 
as well as their respective 95% confidence intervals from Jo = 1 to Jo = 14. We propose to 
choose Jo = 9 for Thul300 and Jo = 10 for Satl930 since from these values of Jo the successive 
confidence intervals are such that the smallest one is included in the largest one (for the robust 
estimators). Note that [23] chose the same values of Jq using another methodology. For these 
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values of Jo we obtain d n: ch = 0.43 (with 95% confidence interval [0.412, 0.443]) , d n ,CR = 0.37 
(with 95% confidence interval [0.358, 0.385]) and <i nj MAD = 0.38 with (95% confidence interval 
[0.362, 0.397]) for Thul300 and (f n>CL = 0.36 (with 95% confidence interval [0.345, 0.374]), 
d ri:CR = d nj MAD = 0.38 (with 95% confidence intervals [0.361, 0.398] for CR and [0.357, 0.402] 
for MAD) for Satl930. These values are similar to the one found by [23]. 





FIGURE 5. Confidence intervals of the estimates d n ,Ch (red), d n ,CR (green) and d n ,MAD 
(blue) on the data Thul300 (left) and Satl930 (right) for Jo = 1, ... , 14 and Jo + £ = 19. 



With this choice of Jo for Thul300, we observe a significant difference between the classical 
estimator and the robust estimators. Thus to better understand the influence of outliers on 
the estimated memory parameter a new dataset with artificial outliers was generated. The 
Thul300 time series shows two spikes shooting down. Especially, the first downward spike 
hits zero. [22] have shown that this dropout lasted 8 seconds. Outliers are introduced by 
dividing by 6 the 8000 observations in this period. The new memory parameter estimators are 
d n ,CL = 0.445, d n ,CR = 0.375 and (in, mad = 0.377. As for the Nile River data, the classical 
estimator was affected while the robust estimators remain stable. 

6. Proofs 

Theorem 4 is an extension of [3, Theorem 4] to arrays of stationary Gaussian processes 
in the unidimensional case and Theorem 5 extends the result of [11] to arrays of stationary 
Gaussian processes. These two theorems are useful for the proof of Proposition 1. 

Theorem 4. Let {Xji, j > l,i > 0} be an array of standard stationary Gaussian processes 
such that for a fixed j > 1, (Xj j)j>o has a spectral density fj and an autocorrelation function 
Pj defined by Pj(k) = K(Xj oXj k), for all k > 0. Assume also that there exists a non increasing 
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sequence {uj}j>i such that for all j > 1 

sup \fj(X) -5oo(A)| < uj , 

AG(-7T,7r) 



(34) 



where is a 2ti -periodic function which is bounded on (— 7r, vr) and continuous at the origin. 
Let h be a function on M with Hermite rank r > 1. We assume that h is either bounded or 
is a finite linear combination of Hermite polynomials. Let {nj}j>\, be a sequence of integers 
such that nj tends to infinity as j tends to infinity. Then, 



rij 

^Y^hiX^^M^a 2 ) , asj 



oo 



(35) 



8=1 



where 



n->oo \^nU*-^ / l\ 

V J i=i e> T 

In the previous equality, q = E,[h(X)Hi(X)], where Hi is the t-th Hermite polynomial and X 
is a standard Gaussian random variable. 



Proof of Theorem 4- Let us first prove that 



A/"(0, 1) , as n — > oo . 



(36) 



'Var (E?i i E,> T ff^fe 
Using Mehler's formula, see Eq. (2.1) of [7], we have 



MEE>(^ 



, i=l 1>t 



E E f^mh(xj,iM 2 (xj,i 2 )} 



11,12=1 h,h>r 



E * 



«2 - »1, 



»i,i2=l 



In order to prove (36), it is enough to prove that for p > 1, 

2 P +r 



E 



(E£l£i>rtffl(*i,i 



0. as n — )• oo and 



(37) 
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For all m G N* , 

rij 



E 



i=l «>r 



E E 



»1- • • • 'm- 



E[^ 1 (X J - il ),...,^ m (X^J] . 



l<il,...,i m <nj /i,...,Z m >r 

1) We start with the case where m = 2p + 1. 

a) Let us first assume that \{ix, ■ ■ ■ , i2 P +i}\ = 2p +1 and that 

Vi, PjW <P* < V(2p) • 



(39) 



By [32, Lemma 3.2 P. 210], E [H h (X jM ),..., H lm (X jA J] is zero if h + • • • + l m is odd. 
Otherwise it is bounded by a constant times a sum of products of (l\ + • • • + i m )/2 corre- 
lations. Bounding, in each product, all of them but p + 1, by p* < l/(2p), we get that 
E [ii^ (Xj^ ),..., H\ 2+1 (Xj t i 2p+1 )~^ is bounded by a finite number of terms of the following 
form 

Z 1 H l~^2p+ 1 / i \ 

(/»*) 3 (P+ Vi(»2-*i)Pi(*4-*3) • • • Pj(i2p-i2p-i)pj{i2p+i-i2 P ) |E {H h {X) . . . H hp+1 (X)) \ 

where X is a standard Gaussian random variable. Note also that the hypercontractivity [32, 
Lemma 3.1 P. 2 10] yields 

|E [H h (X)... H l2p+1 (X)} | < (2p) h+ - + 2 hp+1 V?i!.../ 2p+ i! • 

Thus, using the Cauchy-Schwarz inequality and that p* < there exists a positive constant 
C such that 



E 



Zl,...,Z2p+l>T 



\ c h ■ ■ ■ c h P +i I / * 
III . . . l2p+l ] - 



:- 1± ^-^\K(H h (X)...H hp+1 (X)) 



2p+l 



\ c h I • • • R 2p+ i I ro^„n ' 1+,,, t' 2p+1 -(H 



h,..,h P+1 >r V / ^ L - J 2 P +1 



:(2 PP *)- 



1 )<(2pp*)- 1 (£M [( 2p p *)] 



I P_ 

2 2p+l 



2p+l 
2 



2p+l 
2 



1>1 



^ C El E(W 



2p 
2p+l 



< CO . 



Z>r 



To conclude the proof of (37), it remains to prove that 

Y, Pj{h - h)pj(U - h) ■ ■ ■ Pj{hp - i2 P -i)pj(i2p+i - i2p) 

l<ii,...,i2 p +i<rij 
\{i 1 ,...,i 2p +i}\=2p+l 



0, as rij —7- oo 



£ p\{i2-h) 



H,«2 = l 



(40) 
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Let us first study the numerator in the l.h.s of (40). 



E Pj(i2 - k)Pj(k - h) ■ --Pjihp - i2p-i)Pj(i2 P +i - i2 P ) 



l<ii,...i2 p +i<rij 
\{h,...,i 2p +i}\=2p+l 



= ( E Pj( i 2-h)Y Y Pj{i2p-i2p-l)pj{i2p+l-i2p) 

l<il^i2<rij I<i2p-l,i2p,*2p+i<«j 

|{«2p-l,*2p,*2p+l}|=3 

= ( E Pi( i 2-n)) P ~ 1 E ( E Pj(i2 P +i - kp)) 2 ■ 
To prove (40), we start by proving that 

rij 

Y{ Y Pj(r-s)) =0( nj ). (41) 

r=l l<s<rij 

Using the notation D nj (A) = X^r=l e ' Ar > we 

E( Y P#-*)) 2 = E(f ^ E ^WdA; 

r=l J ~ n Ks<n, 



Ay (A - A') J D„ J (A)A lj (A')/i(A)/ i (A , )dAdA / . 

7T — 7T 



Using (34), the boundedness of g^ and that Uj is bounded, there exists a positive constant C 
such that 

|/i(A)/,(A')| < \fj(X) - 9oo(X)\\fj(X') - ffoo(A')| + booCAOH/iCA) - 5oo (A)| 

+ Iffoo (A) 1 1 /j (A') - goo (A') | + boo (A) | | 5 oc (A') I < C . 

Then, using that there exists a positive constant c such that |_D n . (A)| < aij/(l + n j | A | ) , for 
all A in [— ir, ir], 



rij 



(42) 



The result (41) thus follows from the convergence of the integral in (42) which is proved in 
Lemma 9. Let us now prove that 

— E Pj( r - s ) ~> 5oo(0) , as n ->• oo . (43) 
rij 
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rij 

Using that Fj denned by Fj (A) = ^vm-,-)" 1 ] £ e iAr | 2 , for all A in [-vr, n] satisfies Fj(A)dA = 

r=l 

1, we obtain 

-( E /^- s ))-2oc(0)= r (f,(X)-g OD (X))F j (X)d\+ r (g OD (X)-g oo (0))F,(X)d\. 

(44) 



fc 7 



Using that i ? J (A)dA = 1 and (34), the first term in the r.h.s of (44) tends to zero as n 
tends to infinity. The second term in the r.h.s of (44) can be upper bounded as follows. For 
< rj < 7r, 



< 



\oc(A)- 5 oo(0)|F i (A)dA 



(5oc(A)-5oo(0))F J (A)dA 

+ / |3oo(A)-5oc(0)| J P J (A)dA+ ! | 5oo (A)- 5 oo(0)|i^(A)dA. (45) 

J— T] J T] 

Since there exists a positive constant C such that Fj(X) < C/(?ij|A| 2 ), for all A in [— 7r,7r], the 
first and last terms in the r.h.s of (45) are bounded by Ctt / '(njf] 2 ). The continuity of at 
and the fact that Fj(X)dX < Fj(X)dX = 1 ensure that the second term in the r.h.s of 
(45) tends to zero as n tends to infinity. This concludes the proof of (43). 

Using the same arguments as those used to prove (43) and the fact that pj is the autocor- 
relation associated to fj 1 which is the Z-th self-convolution of fj, we get that 

nj 

~ E P l A r ~ s ) asn^oo. (46) 



n, 

J r,s=l 



Let us now prove that the denominator in (40) is 0(n • 2 ) as n — > oo. We aim at applying 
Lemma 12 with f n , g n , f and g defined hereafter. 

2 I I 

fnM l ) = jM\s\<n J }( 1 -^)^- 

Observe that \f n (s,l)\ < g n \s,l) where 



2 I I 

9nMl) = |l {|s|<n . } (l-l£l)p2 (s) 



Using (34) and the fact that the spectral density associated to p l - is /* , we get, as n — > oo, 

f nj (s,l) -> f(s,l) = | £ g£(X)e iXs dX and g nj (s,l) -> g(s,l) = | £ <4 2 (A)e iAs dA . 
Using [20, Lemma 1], we get 

1>T SGZ 1>T 
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Then, Lemma 12 yields 

■Varf 

rwoo jia 



J 1=1 i>r J 1>t 11,12=1 i>r 

Hence we get (40) by noticing that the numerator in (40) is (9(np. 

If Condition (39) is not satisfied then let k$ be such that Pj(k) < p* < l/(2p), for all k > k^. 
In the case where h is a linear combination of L Hermite polynomials, the same arguments as 
those used previously are valid with p* = 1. In the case where h is bounded, there exists a 
positive constant C such that 



E 



\ 2p+l-| 

K x 3,i)) <C £ E[\h\(X j , ll )...\h\(X j , lq )] , (47) 



n j \ 2p+l 

E 



where ii, . . . , i q are such that \ik — i[\ > ko, for all k, I in {1, . . . , q} with q < 2p + 1. By 
expanding onto the basis of Hermite polynomials, we can conclude with the same arguments 
as those used when Condition (39) is valid. 

b) Let us now assume that ■ ■ ■ , h P +i}\ = r < 2p. In the case where h is bounded, the 
inequality (47) is valid with q < r which gives that the numerator of (37) is 0(nj T ^). In the 

case where h is a linear combination of L Hermite polynomials, we use the same arguments 

lr 
1 j 



as those used in a) with p* = 1 which implies that the numerator of (37) is 0(rv? 



2) Let us now study the case where m is even that is m = 2p with p > 1. 



E 



1 = 1 1>T 



E E ^r^ E [^(^u)...^(^ 

l<ii,...,i'2 P <rij li,...,l2 p >T 1 2p 



2p ■ 



By [26, Formula (33), P.69], we have 



E[H h (X j!h )...H l2p (X j!i2p )]=h\...l 2p \ £ £p (49) 

{Zl,...,i2j>} 

where it is understood that pj = p'j^i.Q ~ k), v\ = v q,k^ an d EE^ i 2 } 

l<q<k<2p l<q<k<2p ' l-i 2p 

indicates that we are to sum over all symmetric matrices v with nonnegative integer entries, 
ua = and the row sums equal to h, . . . , hp- 

We shall prove that among all the terms in the r.h.s of (49), the leading ones correspond 
to the case where we have p pairs of equal indices in the set {li, . . . , hp}, that is, for instance, 
h = h, h = k, ■ ■ ■ , hp-i = hp and z^i ;2 = h, ^3,4 = h,—,^2p-i,2p = h P -i the others Uij being 
equal to zero. This gives 

/, ,N 2 n \\2 Pi( i2 ~ j i)'Vj(»4 ~ h) U ■ --Pjihp ~ i2 P -i) hp 

[h ] -) ■ ■ ■ {hp'-) r-: -r- ; • 

h- ■ ■ ■ hp- 
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The corresponding term in (48) is given by 



clc? ...cl 



;vrV 2 w i2 - ^ l2p ^ H - ^ k ■ ■ ■ rib* - i ^ hp 



l<ii,...,i2 P <nj l2,l4,--;hp>r 



9 J 



= [Y, c i( £ Afc-v) 

1>t «i,i2=l 

which corresponds to the denominator in the l.h.s of (38). Since there exists exactly (2p)\/(2 p p\) 
possibilities to have pairs of equal indices among 2p indices we obtain (38) if we prove that 
the other terms can be neglected. 
Let us first consider the case where 

Vi, Pj (i) < p* < —L- (50) 

and | {ii,..., t2p}| = 2p. By [32, Lemma 3.2 P. 210], E [H h (X jM ),..., H lm (X j>im )] is zero 
if l\ + • • • + l m is odd. Otherwise it is bounded by a constant times a sum of products 
of (Ji + • • • + l m )/2 correlations. Bounding, in each product, all of them but p + 1, by 
p* < l/(2p — 1), we get that E [Hi^Xj^), . . . , H\ 2 (Xj^)] is bounded by a finite number of 
terms of the following form 

(p*) !l± "^ 2E - (p+1 W*2 - h)Pj{H -is)--- Pj{i* P ~ i2 P -i)Pj(i2 P ~ h) |E (H h (X) . . . H hp (X)) \ . 

where X is a standard Gaussian random variable. Using the same arguments as in the case 
where m was odd, we have 

E l ?V""; Q2 ;' ^) !i± ^ £ " (p+1) ! E ( H h w • • • H ^M)) \<oo. 

To have the result (38), it remains to show that 

Yll<ii,-,i2p<nj Pj(h - ii)Pj(,U ~h) ■■■ Pj(h P - i2p-i)pj(fap ~ h) 

l{n '-' l2p}l=2p > , as n -)> oo . (51) 



The numerator of (51) can be rewritten as 

Pji}2 - k)Pj(k - h) ■ ■ ■ Pj(hp - i2p-i)Pj(i2 P - h) 



l<h,...,i2p<nj 
\{il,...,i 2p }\=2p 



( E Pji^-^y E Pj(h-ii)pj(i2p-i2p-i)pj(i2 P -h) 



l<h^ii<nj L I<ii,i2,i2p-i,i2p<rij 

{«l,«2,«2p-l,*2p}|=4 
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Using (43), we have (Xa< i3 ^j 4 < nj Pj(^ — *3) ) = 0(n p - ). Let us now prove that 

Pi(i2 - h)Pj(h - u)pj{k - h) = 0(rij) . (52) 

l<il,j 2 ,«3i*4<«j 

Using the notation D nj (\) = Yl r =l e ' Ar ' 

^2 Pj{i2-h)pj(h-u)pj{h-ix) 

l<il, 22,13, 14<nj 

^ ( r e iA fe-^)/ J (A)dA) ( r j^-^f^w) ( r e^-^/,^) 

l<ii,<a,*s,i4<r»j J ~* J ~ w J w 

~ £(£j D "M\ D nA x + 0\fjW\) /;(£)d£. 

Using (34) and that is bounded, (52) will follow if we prove that f_ (f_ \D nj (X)\\D nj (X + 
£)|dA) 2 d£ = 0(rij) . Since there exists a positive constant c such that |D n -(A)| < crij/(l + 
nj|A|), for all A in [— ir,ir], 

\ 2 poo / poo -j 1 

|D nj (A)||D nj (A+OI/,(A)dA f 3 (Z)6{<c 4 n 3 / d// 



The result (52) thus follows from the convergence of the last integral in (53) which is proved 
in Lemma 10. Hence we get (51) since the numerator of the l.h.s of (51) is 0{n p ~ l ) and 
the denominator is 0{n p -) by the same arguments as those used to find the order of the 
denominator of (40). If Condition (50) is not satisfied or if . . . , ?2p}| < 2p, we can use 
similar arguments as those used in l)a) and l)b) to conclude the proof. □ 

Theorem 5. Let {Xji, j > l,i > 0} be an array of standard stationary Gaussian processes 
such that for a fixed j > 1, (-Xj,i)i>o has a spectral density fj and an autocorrelation function 
Pj defined by pj(k) = K(XjfiXj^), for all k > 0. Let Fj be the c.d.f of Xj t i and F n . the 
empirical c.d.f computed from Xj t i, . . . ,Xj jnj . If Condition (34) holds, 

^n-(F nj -F 3 )^W m D([-oo, oo]) , (54) 



where W is a Gaussian process and D{\— oo, oo]) denotes the Skorokhod space on [—00,00]. 
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Proof of Theorem 5. Let Sj(x) = n- Yli=i \^-{Xj i<x} ~ Fj( x i) ; f° r an x m We shall 
first prove that for x±, . . . , xq and ai, . . . , ag in R 

^a^-^) A AT ^0,^|^(0) j , asrwoo, (55) 

where q is the Z-th Hermite coefficient of the function h defined by 

Q 

^(O = ^a 9 (l{.<xJ-E(l { .< x . g} )) . 

9=1 

Thus, Y^q=i a qSj(x q ) = n- Yli=i h(Xj : i), where h is bounded and of Hermite rank r > 1 
since for all t in R, K(Xlx<t) = f^ x ^x<t l f( x )d- x = J t - 00 {~ l -P{ x ))' & x = ~ L f{t) 0) an d the 
CLT (55) follows from Theorem 4. 

Let us now prove that there exists a positive constant C and (3 > 1 such that for all 
r < s < t, 

E (\Sjis) - S,(r)| 2 |5;(*) " ^(*)| 2 ) < C|t - rf . (56) 
The convergence (54) then follows from (55), (56) and [6, Theorem 13.5]. Note that 

E(|5,( S )-^(r)| 2 |5,(t)-^( S )| 2 ) 

rij rij 

= E E ((/- - fr)( X i,i)(fs ~ fr)(Xj,Mt ~ f.)(X jtl ){f t - fs){X hV ) , 

U i i,i'=ll,l'=l 

where ft(X) = t{x<t} ~ ^(l-TX"<t})- By developing each difference of functions in Hermite 
polynomials , we get 



e (\s 3 ( S ) - swm) - s^i 2 ) =^iziz E 

3 i,i' = ll,l' = l P1 ,...,p 4 >l 

c pi{fs ~ fr)c p2 (f s — fr)c p3 (ft — fs)c p4 (ft — fs) w / tt fv n „ tv \tt \tt Iv \\ 
j j ^ \ H Pi\Xj,i)^p 2 {X j ^)H p3 {X^i)H p4 {X j)V )) . 

pr- - - -P4 ] - 

Using the same arguments as in the case where m is even in the proof of Theorem 4, we obtain 
E(\S j (s)-S s W 2 W)-S J (s)\*) = \ £ jr fv^-f^Ws-U) 



3 pi,p2>ll,l',l,l'=l 



p P^i! _ i)pPf(l> _ /) + C PiU fs)c pi (fs f'Mft fsMfs !A pPf(l _ i)pf(l> _ 

Pl\p2- J j J 
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Let || • H2 = (E(-) 2 ) 1 / 2 and (/, g) = M[f(X)g(X)], where X is a standard Gaussian random 
variable. Since, by (46), ^ nr =1 f^ 1 (I — %)(?^{V — i') = 0(re 2 ), we get with the Cauchy- 
Schwarz inequality that there exists a positive constant C such that 

E (|,S,(„) - S;(r)| 2 |S 3 -(t) - S,( a )\ 2 ) 



PliP2>l 



c li ~ f s "> C P2 (/« ~ fr) , Cpi (/t - /s)Cpi (/s - /r)c p2 (/t - / s )c p2 (/ s - /, 



"P2 



+ 



< C (\\f t - fs\\l \\fs - fr\\l + \(ft- fs, fs ~ fr)?) < C Hf t - f s \\l \\f s - f 



2 

r\\2 



Note that lift- fall! < 2(||l { x<t} 

|2 



1 {^<4ll2 + ll E ( 1 {^<4)- E ( 1 {^})|l2)- SillCeS ^ *> 

l{x<t} ~~ l{x<s}||2 = <&(t) — $(s) < C\t— s |, where $ denotes the c.d.f of a standard Gaussian 
random variable. Moreover, ||E(l{x<s}) ~~ ^(l{X<t})||2 — s | 2 > which concludes the proof 
of (56). □ 

Proof of Proposition 1. We first prove (21) for * = CL. 



1 nj ~ L 2cr 2 nj ~ 1 1 /tf/ 2 



(7,- 



i=0 



J i=o J 



Let us now prove (21) for * = MAD. Let us denote by F n . the empirical c.d.f of Wj t o :n .-i and 
by Fj the c.d.f of W jj0 - Note that 

^MADj = rn($)T (F n .) , 

where T = T 2 o Ti with Ti : F H> {r H> / R l { | x .|< r} dF(2;)} and T 2 : U ^ Z7 _1 (l/2). To 
prove (21), we start by proving that ^/nj(F nj — Fj) converges in distribution in the space 
of cadlag functions equipped with the topology of uniform convergence. This convergence 
follows by applying Theorem 5 to Xji = Wj t i/(Tj which is an array of zero mean stationary 
Gaussian processes by [20, Corollary 1]. The spectral density of (Xji)i>o is given by 
fj(X) = Dj j o(A; f)/(Jj where D J: o(-;/) is the within scale spectral density of the process 
{Wj,k}k>o defined in (12) and cr 2 is the wavelet spectrum defined in (14). Here, <7oo(A) = 
D OOj0 (A;(i)/K(d), with D OOi0 (-;d) defined in (13) and K(d) = \^ 2d m)\ 2 ^ since by 
[20, (26) and (29) in Theorem 1] 

D i)0 (A;/) D oo , (A;d) 



f*(0)K(d)2 2d i K(d) 



/*(0)K(d)2 2 *' 



< CLK(d)- 1 2-^ j 
< CL2~ l3j ->-0, asn^oo. 



, as n — > 00 , 
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Note also that, by [20, Theorem 1], <?oo(A) is a continuous and 2-7r-periodic function on {—it, it). 
Moreover, g oc (X) is bounded on (— tt,tt) by Lemma 11 and 

where Cx and C2 are positive constants. The asymptotic expansion (21) for ctmaDj can be 
deduced from the functional Delta method stated e.g in [33, Theorem 20.8] and the classical 
Delta Method stated e.g in [33, Theorem 3.1]. To show this, we have to prove that To = T10T2 
is Hadamard differentiable and that the corresponding Hadamard differential is defined and 
continuous on the whole space of cadlag functions. We prove first the Hadamard differentia- 
bility of the functional T±. Let (gt) be a sequence of cadlag functions with bounded variations 
such that \\gt — g\\oo — > 0, as t — > 0, where g is a cadlag function. For any non negative r, we 
consider 

T^+tg^rj-T^jr] _ (Fj + tg t )(r) — (Fj + tg t )(-r) — Fj(r) + Fj{—r) 

t t 

tgt(r) - tgt(-r) 
= 7 = 9t{r) - gt(-r) -»• g{r) - g(-r), 

since \\gt — g\\oo — > 0, as t — > 0. The Hadamard differential of T\ at g is given by : 

(DTi(Fj).g)(r) = g(r) — g(—r). 

By [33, Lemma 21.3], T2 is Hadamard differentiable. Finally, using the Chain rule [33, The- 
orem 20.9], we obtain the Hadamard differentiability of Tq with the following Hadamard 
differential : 

nT(F) (DT^.g^TojFj)) g(T (Fj)) - gj-TojFj)) 

0Kl) - 9 - (T^FjmUFj)] ~ (Ti(Fj))'[T (Fj)] ' 

In view of the last expression, DTq(Fj) is a continuous function of g and is defined on the 
whole space of cadlag functions. Thus by [33, Theorem 20.8], we obtain : 

m($)^n] (T (F nj ) - T (F,)) = m{$>)DT G (F 3 ) {^](F nj - Fj)} + o P (l), 

where m(<I>) is the constant defined in (19). Since To(Fj) = <Tj/m($) and (Ti(Fj))'(r) = 
2aj~ 1 ip(r /dj), where 99 is the p.d.f of a standard Gaussian random variable, we get 

<Jnj{duADj-<rj) = -7= IF ( — ,MAD,$ ) +o P (l) 

and the expansion (21) for * = MAD follows from the classical Delta method applied with 
f(x) = x 2 . We end the proof of Proposition 1 by proving the asymptotic expansion (21) for 
* = CR. We use the same arguments as those used previously. In this case the Hadamard 
differentiability comes from [16, Lemma 1]. □ 
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The following theorem is an extension of [3, Theorem 4] to arrays of stationary Gaussian 
processes in the multidimensional case. 

Theorem 6. Let Xj^ = {Xj*J , . . . , Xj} } be an array of standard stationary Gaussian pro- 
cesses such that for in {0, . . . ,d}, the vector {XjJ,Xjfy has a cross-spectral density 

(i i') (i i') Hi 1 )/ \ i (i) (i 1 ) \ 

fj and a cross-correlation function pj defined by pj {kj = ~E(Xj:Xj i+k ), for all 

k > 0. Assume also that there exists a non increasing sequence {uj}j>i such that uj tends to 

zero as J tends to infinity and for all J > 1, 

sup \ri j ' j 'H\)-g&'\\)\<uj, (57) 

Ae(-7r,7r) 

where ^ is a In -periodic function which is bounded on (—it, 7r) and continuous at the origin. 
Let h be a function on R with Hermite rank r > 1 which is either bounded or is a finite linear 
combination of Hermite polynomials. Let (3 = {/3q, . . . , /3d} in M. d+1 and T~L : — > M. the 
real valued function defined by %(x) = ^j=Q&jh(xj). Let {nj}j>i be a sequence of integers 
such that nj tends to infinity as J tends to infinity. Then 



where 




Ln the previous equality, C£ = E,[h(X)Hi(X)], where Hg is the i-th Hermite polynomial and X 
is a standard Gaussian random variable. 

The proof of Theorem 6 follows the same lines as the one of Theorem 4 and is thus omitted. 

Proof of Theorem 2. Without loss of generality, we set /*(0) = 1. In order to prove (25), let 
us first prove that for a = (ao, • • • , ott) where the ctj's are in M, 
e 

V / ^o 2 - 2J ° d 5> j (al Jo+j (Wj 0+j ,o:n Jo+j -l) -<j 0+i ) AAT(0,a T U»(rf)a) . (59) 
j=0 

By Proposition 1, 

I 

3=0 

^^Vi t) EV^,*) t ,d,. (60) 
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Thus, proving (59) amounts to proving that 

nj +j-l 

y 



since ct 2 o+j . y/n2- J °2~ 2Jod /n Jo+j ~ 2 2 ^/ 2+ ^K(d)/*(0)/ 1 /n^ , as n tends to infinity, by 
[20, (29) in Theorem 1]. Note that 



nj 0+^ 1 /rrr \ ^0+^2^-1 , rrr 



nj n +i— 1 



+ J y wf^M.,^ 



Using the notation: /3, = 2ay2 2d:? ^/ 2+J 'K(d)/*(0) and that IF is either bounded or equal to 
1 fl it? I W Jo+j,i 



Eft E IF (^^J 



i n Jo+e -l 2 e-j_ 1 



y F(Y JoAi ,*) + op(l), 



where 



and 



F(Y JoAu *)=± Pj 2 y 1 IF (eWiI^, *, $ 



j=0 u=0 



Wj 0+ £ ti Wj 0+ g-i t 2i Wj 0+ e-i t 2i+l W Jo+j2 e-H 



, . . . , 



^J +i,2 f -ij+2 f -i-l W J ,2 e i W J ,2 e i+2 e -l 



is a 2^ +1 — 1 stationary Gaussian vector. By Lemma 7, F is of Hermite rank larger than 2. 
Hence, from Theorem 6 applied to %{■) = F(-), Xjj = Yj ,e,i an d fo(-) = IF(-), we get 

nj +*-l 

-— = ^ F(y 7o/ „*) AAA(0,5 2 ), (62) 



28 



O. KOUAMO, C. LEVY-LEDUC, AND E. MOULINES 



where a'i = lim^ nJ^Var [ £> +< * F (Xj ,i,U *)) ■ By [20, (26) and (29)] and by using 
the same arguments as those used in the proof of Proposition 1, Condition (57) of Theorem 6 
holds with ff j '\\) = H /.,,,., y{\:f) : a. h ,. l a. h .y and gg^ = T}^_ .,(A; d)/K(d), where 
< r < 2 J ~ J — 1 and Dj 0+ jj_j/(-; /) is the cross-spectral density of the stationary between 
scale process defined in (12). Lemma 11 and [20, Theorem 1] ensure that T^j_.,{-\d) is a 
bounded, continuous and 27r-periodic function. 

By using Mehler's formula [7, Eq. (2.1)] and the expansion of IF onto the Hermite polymials 
basis given by: IF(x,*,$) = £ P >2 c pC IF *)- H p( a; )M where c p (IF*) = E[IF(X, *, $)H p (X)], 
H p being the pth. Hermite polynomial, we get 



1 , n ^~ L N 

— Varf Y, F (YJo^*)) 



1 



nj0+£ .i,>=l 



E E E E E 



i,i'=0 v=0 v'=0 



IF 



W 



J 0+j ,2'-H +V \ w 



a 



Jo+j 



w 



Jo +j'^-j'i'+v' 
a Jo+j' 



*, $ 



— E Wr E E e[if(1^,*,*)if(5^ j *,*)1 +0 (i) 



i=0 i'=o 

n J +3- ln J +3 , ~ 1 



a Jo+j' 



nj 0+ i .-7 =1 



E«" E E E 



+ o(l). (63) 



i=0 i'=0 p>2 

Without loss of generality, we shall assume in the sequel that j > f. (63) can be rewritten as 
follows by using that i' = 2 j ~i q + r, where g 6 N and r € {0,1,... , 2^ - 1} and Eq. (18) 
in [20] 



^ E m E E E E» 



i=0 <j=0 r=0 p>2 
2i-i'-l 



P! 



^E«' E E E^^ 



|T|<nj Q+J r=0 p>2 



P 



E 



n Jo+j / l cr Jo+i 



Wj 0+ j,Q W J +j>,2i-i'T+r 



°~Jo+j' 



+o(l) 



n Jo+j 



•23-3 _1 



E m E E E 



<£( IF *) ^ |r 



pi 



1- 



VJa+jVJo+j' 



dA +o(l) 



|r|<nj 0+i r=0 P>2 

where Dj 0+JiJ _j'(-; /) is the cross-spectral density of the stationary between scale process 
defined in (12). We aim at applying Lemma 12 with f n , g n , f and g defined hereafter. 



fnj +j ( T i P) 



P i 



E l{M<n W (l-^ 



E 



w Jo+j , w Jo+j 



',23-3 T+r 
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Observe that \f njo+j \ < 9n Jn+i , where 



''Jq+3 ' 



fj, 



Jo+j 



pi 



H\r\<nj 0+j }(l 



Using [20, (26) and (29) in Theorem 1] we get that 

'>/,.,., ,'iA:/i _ 2 d ^') 



E 



a Jo+j 



+j,0 yy Jo+j',2i-i T+r 



lim 

n ^°° VJ +jVJ +j> 



K(d) 

This implies that lim^^ f nj Ar,p) = /(r,p) where 



Doo,j-j'(A; d) 



f(r,p) 



23-3 _1 

E 

r=0 



gg^2 r d« 

K(d) 



,(A;d)e iAr dA 



00,3-3 



Futhermore, Umn^^ g njo+ .(r,p) = g{r,p) where 



9(t,p) 



c 2 p (lF*) 2 M (i-i') 



p! K(d) 2 



Doo,j-j'(A;d)e iAr dA 



and |x|| = Y7k=i xf, for x = (x\, . . . , x r ) G R r . Using (63)-(65) in [20] we get 



EE 



9i 



Ja+i 



p>2 T ez 

Then, with Lemma 12, we obtain 



23-3 _1 



? . 2 = E C ' { ™mui l T E 4c,« j ,2*<^)2*'P-pH/ E £ 



p!K(d)P- 2 
P >2 ^ v y jii/=0 



Dg i _.,(A;d)e^dA 

□ 



c 2 (IF*)\ 2 2d ( j ~ j,S) 
S P p ! J K(rf) 2 27F / |Doo,i-i'(A;d)||dA 



7. Technical Lemmas 

Lemma 7. Lei X 5e a standard Gaussian random variable. The influence functions IF defined 
in Proposition 1 have the following properties 



E[IF(X,*,c&)] = 0, 
E[XIF(X, *, $)] = , 
E[X 2 IF(X ,*,$)] + . 



(64) 
(65) 
(66) 



Proof of Lemma 7. We only have to prove the result for * = MAD since the result for * = CR 
follows from [16, Lemma 12]. (64) comes from E(l{x<i/m(*)}) = ■^(l{X<$- 1 (3/4)}) = 3/4 and 
^(1{X<— l/W$)}) = l/^) where X is a standard Gaussian random variable. (65) follows from 
I R xt {x ^-i {3/4)} <f(x)dx - J H j;l {:r <_ <I> -i ( 3 /4)} ^(a;)d2; = -p($ _1 (3/4)) + ^(-$ _1 (3/4)) = 0, 
where <p is the p.d.f. of a standard Gaussian random variable and the fact that E(X) = 0. Let 
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us now compute E[X 2 IF(A, MAD, $)]. Integrating by parts, we get J R x 2 l^ x< ^-i^/^jip(x)dx— 
3/4-/ R x 2 l {x <_ $ - 1(3/4)} ^(x)dx + l/4 = -2<p i*-^/*)) • Thus, E[X 2 IF(X, MAD, $)] = 2 ± 
0, which concludes the proof. □ 

Lemma 8. Let X be a standard Gaussian random variable. The influence functions IF defined 
in Lemma 1 have the following properties 

E[IF 2 (X, MAD, *)] = ^^(^ = 1-3601 , (67) 

E[IF 2 (X, CR, $)] » 0.6077 . (68) 

Proof of Lemma 8. Eq (68) comes from [27]. Since , 

E[IF 2 (A,MAD,$)] = i -^g|^— Var (l {W <*-i(s/4)}) , 
where l{|x|<<l>- 1 (3/4)} is a Bernoulli random variable with parameter 1/2, (67) follows. □ 
Lemma 9. 

Proof of Lemma 9. Let us set I = p(/j,, /J,')d/id(i' , with 

( 1 1 1 

^'^-l + lM-M'll + Hl + lM'l ' 

Note that I = /i + J 2 + ^3 + h, where h = J °° / °° //)d/xd//, J 2 = Jo°° /"^ p(jJ>, /t/)d/xd//, 

I3 = p(/J,, fj,')dfj,dfi' and J 4 = $_ 00 p(,l JI >) //)d/id//. It is easy to see that Ii = I4 and 

/ 2 = ^3- Let us now compute Ji. Using partial fraction decomposition, 

r°° 1 / r°° 1 1 \ [°° 1 / 1 1 \ 

r iog(i+^o , , , 9 r iog(i+^) i , 

= /o 7uT77 d ^ +2 /o (2 + ,o(i + .') dM<00 ' 

since in the neighborhood of 0, log(l + //)/{//(l +aO) ~ V^ + Z-Oj log(l + //)/{(2 + //)(1 + 
//')} ~ {-1/(1 + //) +2/(2 + //)} and in the neighborhood of 00, log(l + //)/{//(l + //)} and 
log(l + //)/{(2 + //)(l + //)} ~ log(//)/ // 2 . Let us now compute I 2 . Using the same arguments 
as previously, we get 

h =r-±-(r u„w = r lo f +/ 1W < «, . 



□ 



Lemma 10. 



/.(/. i + J+M-U + IH '") d "' <0 °' 
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Proof of Lemma 10. Let us set / = ( p(/j,, fx')d/j,) d//, with 

1 + |/i + 1 + |/i| 

Note that / < 2(I 1 +/ 2 +J 3 +/ 4 ), where ^ = / °° ( / °° p(/i, //)d/i) 2 d//, h = / °° ( p(m, //)d/i) ; 
4 = J°oo (J^° pip., fjf)dfj) 2 dfj,' and J 4 = /"^ ( /"^(/i, p!)dp) 2 dp!. It is easy to see that 
7l = I 4 and / 2 = ^3- Let us now compute /j. Using partial fraction decomposition, 

1 



l + /i + //l + /i 



1 \ 2 

d/i I d/i' 



o 



00 ' 1 \2 

- log(l + //) d/x' < oo 



since in the neighborhood of 0, [log(l + /i')] 2 //i' 2 ~ 1, and in the neighborhood of oo, [log(l + 
/i')] 2 j fj,' 2 ~ [log(/i')] 2 //i' 2 . Let us now compute 12- Using the same arguments as previously, 
we get that there exists a positive constant C such that 



h < 2 



1 1 \ 2 

d/i) d/i' + 2 



1 + /i - pi 1 + /i / 

'■ oo 



1 



1 - /i + pi 1 + p 



1 \ 2 
d/i ) d/i' 



o 



Mi + iJ_ 



i PC 

dp' + C 

Jo 



mi + d 

(2 + /*') 



d/i' < oo 



□ 



Lemma 11. Lei e u (f) = 2~ u / 2 [l,e~ i2 . . . , e"^ 2 "" 1 ) 2 U ^] T , where (el. For all u > 0, 
each component of the vector 



Doo, u (A;d) = ^|A + 2/7rr 2d e u (A + 2/7r)V(A + 2/7r)V'(2-"(A + 2/7r)) , 
is bounded on (— vr, vr), where ip is defined in (4). 

Proof of Lemma 11. We start with the case where I = 0. Using (5), we obtain that 
2~«/ 2 |A|- M |$(A)|$(2- u A)| = 0(\X\ 2M ~ 2d ), as A -»• hence, (7) ensures that 
2~ M / 2 |A|~ 2d |^(A)||^(2- u A)| = 0(1). Let ei fc) denotes the fe-th component of the vector e u . 
For I ^ 0, (W-2) ensures that for all A in (— vr, vr) there exists a positive constant C such that 
IV'(^)I — C/(l + 1^1)°- Then, there exists a positive constant C" such that 

^ |A + 27rl\~ 2d ^(\ + 2ttZ)^(2- u (A + 2vr/))e^(A) < C ^ |A + 2^r M " 2a . 

lez* lez* 

If A = 0, £ 1/|2vtZ| m+2q < oo by (7). If A / 0, then, since -vr < A < vr, £ 1/|A + 2vrZ| 2d+2a < 

zez* ' zez* 
E l/|vr(2Z-l)| 2d+2Q <ooby (7). 



□ 



Lemma 12. Let f n and g n be two sequences of measurable functions on a measure space 
({l,J-,p) such that for all n \f n \ < g n . Assume thatlimm£g n exists and is equal to g . Assume 

n— >oo 

also that J gdp = liminf J g n dp and lim f n = f. Then J lim inf f n dp = liminf J f n dp. 
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Proof of Lemma 12. Since f n = /+ — / , where /+, f n > 0, we assume in the sequel that f n 
is non negative. By Fatou's Lemma f liminf(g n — f n )dfi < liminf f {g n — f n )dfJ,. Using that 

n— >oo n— >oo 

liminf g n = g and that j gdfi = liminf J g n dfi, we obtain limsup j fndfJ. < Jlimsup/ n d/i. By 

71— >oo n— >oo n— >oo n— >oo 

applying Fatou's Lemma to f n , we obtain f lim inf f n d^ < liminf f f n d/x. Thus, 

n— >oo n— ^oo 

/ fd/j, = / lim inf f n d/j, < lim inf / f n d^i < lim sup / / n d/i < / lim sup/ n d// = / fdfi , 

J J rwoo n-i>oo J J J J 

which concludes the proof. □ 
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